Non-Linear Mixed Models

---

Packages

library(tidyverse)
library(nlme)
library(plotly)

Importing Dataset

setwd("~/Insync/leonardogloria@uenf.br/Google Drive/MENTORIA_R/Monday_meeting")
dadosg <- read.table("data_growth.txt", h = T,na.strings=".") %>% 
  mutate(age=idade,w=peso,animf=factor(animal),class=factor(sexo)) %>% 
  groupedData(w~age|animf,data=.) #Creating dataset indicating animal repeated measures by age

Logistic model, without fixed effect

teta=c();   AIcc=c();

model00=nlme(w~a/(1+b*(exp(-k*age))),fixed=list(a+b+k~1),random=pdDiag(a~1),
             control=nlmeControl(maxIter=100),
             data=dadosg,
             start=list(fixed =c(30.4443882,12.7867281,0.06)),na.action=na.omit)       

(teta[1]=attributes(logLik(model00))$df)
(AIcc[1]=-2*model00$logLik+2*teta[1]+(2*teta[1]*(teta[1]+1))/(nrow(dadosg)-teta[1]-1))

Logistic model with class of fixed effect, in other words, we will fit a parameter for each sex

st01=c(rep(summary(model00)$ coefficients$fixed[[1]],length(levels(dadosg$class))),rep(summary(model00)$ coefficients$fixed[[2]],length(levels(dadosg$class))),rep(summary(model00)$ coefficients$fixed[[3]],length(levels(dadosg$class)))) #starting values based on the model00 output

model01=nlme(w~a/(1+b*(exp(-k*age))),fixed=list(a+b+k~class-1),
             random=pdDiag(a+b+k~1),control=nlmeControl(maxIter=500),
             data=dadosg,start=st01,na.action=na.omit)


(teta[2]=attributes(logLik(model01))$df)
(AIcc[2]=-2*model01$logLik+2*teta[2]+(2*teta[2]*(teta[2]+1))/(nrow(dadosg)-teta[2]-1))

Logistic model with class of fixed effect, in other words, we will fit a parameter for each sex, plus (weights = heterogeneous variance along the time, varPower= power variance function is defined as s2(v) = |v|^(2*t))

parf=c()
for (i in 1:length(summary(model01)$ coefficients$fixed))      {
par1=summary(model01)$ coefficients$fixed[[i]]
parf=rbind(parf,par1)
 } #starting values based on the model01 output, each class has its own start value

model02=nlme(w~a/(1+b*(exp(-k*age))),fixed=list(a+b+k~class-1),
             random=pdDiag(a+b+k~1),control=nlmeControl(maxIter=500),
             data=dadosg,start=parf,na.action=na.omit,weights= varPower())       

intervals(model02)
######################

(teta[3]=attributes(logLik(model02))$df)
(AIcc[3]=-2*model02$logLik+2*teta[3]+(2*teta[3]*(teta[3]+1))/(nrow(dadosg)-teta[3]-1))

Logistic model with class of fixed effect, in other words, we will fit a parameter for each sex, plus (corr = modeling animal repeated measures, CAR1 = continous autoregressive)

parf=c()
for (i in 1:length(summary(model01)$ coefficients$fixed))      {
par1=summary(model01)$ coefficients$fixed[[i]]
parf=rbind(parf,par1)
 } #starting values based on the model01 output, each class has its own start value

model03=nlme(w~a/(1+b*(exp(-k*age))),fixed=list(a+b+k~class-1),
             random=pdDiag(a+b+k~1),control=nlmeControl(maxIter=100),
             data=dadosg,start=parf,na.action=na.omit,corr=corCAR1())
######################

(teta[4]=attributes(logLik(model03))$df)
(AIcc[4]=-2*model03$logLik+2*teta[4]+(2*teta[4]*(teta[4]+1))/(nrow(dadosg)-teta[4]-1))

Full Logistic model with class of fixed effect, heterogeneous variance along the time, repeated measures

parf=c()
for (i in 1:length(summary(model03)$ coefficients$fixed))      {
par1=summary(model03)$ coefficients$fixed[[i]]
parf=rbind(parf,par1)
 }

model04=nlme(w~a/(1+b*(exp(-k*age))),fixed=list(a+b+k~class-1),
             random=pdDiag(a~1),control=nlmeControl(minScale=10**-100,maxIter=500),
             data=dadosg,start=parf,na.action=na.omit,corr=corCAR1(),weights= varPower())       #logistico
######################

(teta[5]=attributes(logLik(model04))$df)
(AIcc[5]=-2*model04$logLik+2*teta[5]+(2*teta[5]*(teta[5]+1))/(nrow(dadosg)-teta[5]-1))

Gompertz model, without fixed effect

model05=nlme(w~a*exp(-b*(exp(-k*age))),fixed=list(a+b+k~1),random=pdDiag(a~1),
             control=nlmeControl(maxIter=100),
             data=dadosg,
             start=list(fixed =c(30.4443882,12.7867281,0.06)),na.action=na.omit)       

(teta[6]=attributes(logLik(model05))$df)
(AIcc[6]=-2*model05$logLik+2*teta[6]+(2*teta[6]*(teta[6]+1))/(nrow(dadosg)-teta[6]-1))

Gompertz model with class of fixed effect, in other words, we will fit a parameter for each sex

st05=c(rep(summary(model05)$ coefficients$fixed[[1]],length(levels(dadosg$class))),rep(summary(model05)$ coefficients$fixed[[2]],length(levels(dadosg$class))),rep(summary(model05)$ coefficients$fixed[[3]],length(levels(dadosg$class)))) #starting values based on the model05 output

model06=nlme(w~a*exp(-b*(exp(-k*age))),fixed=list(a+b+k~class-1),
             random=pdDiag(a+b+k~1),control=nlmeControl(maxIter=100),
             data=dadosg,start=st01,na.action=na.omit)


(teta[7]=attributes(logLik(model06))$df)
(AIcc[7]=-2*model06$logLik+2*teta[7]+(2*teta[7]*(teta[7]+1))/(nrow(dadosg)-teta[7]-1))

Gompertz model with class of fixed effect, in other words, we will fit a parameter for each sex, plus (weights = heterogeneous variance along the time, varPower= power variance function is defined as s2(v) = |v|^(2*t))

parf=c()
for (i in 1:length(summary(model06)$ coefficients$fixed))      {
par1=summary(model06)$ coefficients$fixed[[i]]
parf=rbind(parf,par1)
 } #starting values based on the model01 output, each class has its own start value

model07=nlme(w~a*exp(-b*(exp(-k*age))),fixed=list(a+b+k~class-1),
             random=pdDiag(a+b+k~1),control=nlmeControl(maxIter=100),
             data=dadosg,start=parf,na.action=na.omit,weights= varPower())       

intervals(model07)
######################

(teta[8]=attributes(logLik(model07))$df)
(AIcc[8]=-2*model07$logLik+2*teta[8]+(2*teta[8]*(teta[8]+1))/(nrow(dadosg)-teta[8]-1))

Gompertz model with class of fixed effect, in other words, we will fit a parameter for each sex, plus (corr = modeling animal repeated measures, CAR1 = continous autoregressive)

parf=c()
for (i in 1:length(summary(model06)$ coefficients$fixed))      {
par1=summary(model06)$ coefficients$fixed[[i]]
parf=rbind(parf,par1)
 } #starting values based on the model01 output, each class has its own start value

model08=nlme(w~a*exp(-b*(exp(-k*age))),fixed=list(a+b+k~class-1),
             random=pdDiag(a+b+k~1),control=nlmeControl(maxIter=100),
             data=dadosg,start=parf,na.action=na.omit,corr=corCAR1())
######################

(teta[9]=attributes(logLik(model08))$df)
(AIcc[9]=-2*model08$logLik+2*teta[9]+(2*teta[9]*(teta[9]+1))/(nrow(dadosg)-teta[9]-1))

Full Gompertz model with class of fixed effect, heterogeneous variance along the time, repeated measures

parf=c()
for (i in 1:length(summary(model06)$ coefficients$fixed))      {
par1=summary(model06)$ coefficients$fixed[[i]]
parf=rbind(parf,par1)
 } #starting values based on the model01 output, each class has its own start value

model09=nlme(w~a*exp(-b*(exp(-k*age))),fixed=list(a+b+k~class-1),
             random=pdDiag(a+k~1),control=nlmeControl(minScale=10**-100,maxIter=500),
             data=dadosg,start=parf,na.action=na.omit,corr=corCAR1(),weights= varPower())
######################

(teta[10]=attributes(logLik(model09))$df)
(AIcc[10]=-2*model09$logLik+2*teta[10]+(2*teta[10]*(teta[10]+1))/(nrow(dadosg)-teta[10]-1))

multi-model selection framework (Burnham and Anderson, 2004)

delta=c()
for(i in 1: length(AIcc)){
delta[i]=AIcc[i]-min(AIcc)
}

wpro=c()
for(i in 1:length(AIcc)){
wpro[i]=exp(-delta[i]/2)
}
sum(wpro[1:length(AIcc)])

wprob=c()
for(i in 1: length(AIcc)){
wprob[i]=exp(-delta[i]/2)/sum(wpro[1:length(AIcc)])
}


ER=c()
for(i in 1: length(AIcc)){
ER[i]=max(wprob)/wprob[i]
}

(quadro.akaike=data.frame(teta,AICc=AIcc,delta,wprob,ER))

Extracting the confidence intervals of the parameters after choosing the best model

intervals(model09)

Prediction and Growth rate function

Gompertzhat=function(a,b,k,age){
  y=a*exp(-b*(exp(-k*age)))
  return(y)
}
GRGompertz=function(a,b,k,age){
  y=a * (exp(-b * (exp(-k * age))) * (b * (exp(-k * age) * k)))
  return(y)
}

creating interactive graphics with plotly

setwd("~/Insync/leonardogloria@uenf.br/Google Drive/MENTORIA_R/Monday_meeting")
param=read.table("parameters.txt",h=T)
row.names(param)=c("Males","Females")

time=seq(1,90,by=0.25)
#predictions
predc1=Gompertzhat(a=param[1,1],b=param[1,2],k=param[1,3],age =time)
predc2=Gompertzhat(a=param[2,1],b=param[2,2],k=param[2,3],age =time)

#growth rate
grc1=GRGompertz(a=param[1,1],b=param[1,2],k=param[1,3],age =time)
grc2=GRGompertz(a=param[2,1],b=param[2,2],k=param[2,3],age =time)

#inflaction point (max efficiency)
maxef1=data.frame(grc1)
rownames(maxef1)=time

maxefx1=as.numeric(rownames(maxef1)[which.max(maxef1$grc1)])
maxefy1=maxef1$grc1[which.max(maxef1$grc1)]

maxef2=data.frame(grc2)
rownames(maxef2)=time

maxefx2=as.numeric(rownames(maxef2)[which.max(maxef2$grc2)])
maxefy2=maxef2$grc2[which.max(maxef2$grc2)]

### Multiple Y Axes

ay <- list(
  tickfont = list(color = "black"),
  overlaying = "y",
  side = "right",
  title = "Daily Weight Gain(g)"
)
fig <- plot_ly()
fig <- fig %>% add_lines(x = time, y = predc1, name = "Weight Males")
fig <- fig %>% add_lines(x = time, y = grc1, name = "Growth Rate Males", yaxis = "y2")
fig <- fig %>% add_markers(x = maxefx1, y = maxefy1, name = "Max Growth Rate Males", yaxis = "y2")
fig <- fig %>% add_segments(line=list(type="dash"),x = maxefx1,xend=maxefx1, y = 0,yend =maxefy1, yaxis = "y2",name = "Max Growth Rate Males")

fig <- fig %>% add_lines(x = time, y = predc2, name = "Weight Females")
fig <- fig %>% add_lines(x = time, y = grc2, name = "Growth Rate Females", yaxis = "y2")
fig <- fig %>% add_markers(x = maxefx2, y = maxefy2, name = "Max Growth Rate Females", yaxis = "y2")
fig <- fig %>% add_segments(color = I("gray"),x = maxefx2,xend=maxefx2, y = 0,yend =maxefy2, yaxis = "y2",name = "Max Growth Rate Females")

fig <- fig %>% layout(
  title = "Growth Curve with Growth rate", yaxis2 = ay,
  xaxis = list(title="Age(days)"),
  yaxis = list(title="Weight(g)")
)
fig

Importing Milking Dataset

Importing dataset from txt

setwd("~/Insync/leonardogloria@uenf.br/Google Drive/MENTORIA_R/Monday_meeting")
data_milkg <- read.table("data_milk.txt", h = T) %>% 
  mutate(animal=factor(animal),CG=factor(CG)) %>% 
  groupedData(MY~DIM|animal,data=.) #Creating dataset indicating animal repeated measures by days in milk

WOOD model, without fixed effect

teta=c();   AIcc=c();

model00=nlme(MY~a*(DIM**b)*exp(-k*DIM),fixed=list(a+b+k~1),random=pdDiag(a+b~1),
             control=nlmeControl(maxIter=100),
             data=data_milkg,
             start=list(fixed =c(20,0.1,0.003)),na.action=na.omit)       

(teta[1]=attributes(logLik(model00))$df)

[1] 6

(AIcc[1]=-2*model00$logLik+2*teta[1]+(2*teta[1]*(teta[1]+1))/(nrow(data_milkg)-teta[1]-1))

[1] 13428.63

WOOD model with class of fixed effect, in other words, we will fit a parameter for each sex

st01=c(rep(summary(model00)$ coefficients$fixed[[1]],length(levels(data_milkg$CG))),rep(summary(model00)$ coefficients$fixed[[2]],length(levels(data_milkg$CG))),rep(summary(model00)$ coefficients$fixed[[3]],length(levels(data_milkg$CG)))) #starting values based on the model00 output

model01=nlme(MY~a*(DIM**b)*exp(-k*DIM),fixed=list(a+b+k~CG-1),random=pdDiag(a~1),
             control=nlmeControl(maxIter=100),
             data=data_milkg,start=st01,na.action=na.omit)


(teta[2]=attributes(logLik(model01))$df)

[1] 20

(AIcc[2]=-2*model01$logLik+2*teta[2]+(2*teta[2]*(teta[2]+1))/(nrow(data_milkg)-teta[2]-1))

[1] 13355.63

WOOD model with class of fixed effect, in other words, we will fit a parameter for each sex, plus (weights = heterogeneous variance along the time, varPower= power variance function is defined as s2(v) = |v|^(2*t))

parf=c()
for (i in 1:length(summary(model01)$ coefficients$fixed))      {
par1=summary(model01)$ coefficients$fixed[[i]]
parf=rbind(parf,par1)
 } #starting values based on the model01 output, each class has its own start value

model02=nlme(MY~a*(DIM**b)*exp(-k*DIM),fixed=list(a+b+k~CG-1),random=pdDiag(a~1),
             control=nlmeControl(maxIter=500),
             data=data_milkg,start=parf,na.action=na.omit,weights= varPower())       

intervals(model02)
######################

(teta[3]=attributes(logLik(model02))$df)
(AIcc[3]=-2*model02$logLik+2*teta[3]+(2*teta[3]*(teta[3]+1))/(nrow(data_milkg)-teta[3]-1))

WOOD model with class of fixed effect, in other words, we will fit a parameter for each sex, plus (corr = modeling animal repeated measures, CAR1 = continous autoregressive)

parf=c()
for (i in 1:length(summary(model01)$ coefficients$fixed))      {
par1=summary(model01)$ coefficients$fixed[[i]]
parf=rbind(parf,par1)
 } #starting values based on the model01 output, each class has its own start value

model03=nlme(MY~a*(DIM**b)*exp(-k*DIM),fixed=list(a+b+k~CG-1),random=pdDiag(a~1),
             control=nlmeControl(maxIter=100),
             data=data_milkg,start=parf,na.action=na.omit,corr=corCAR1())
######################

(teta[4]=attributes(logLik(model03))$df)

[1] 21

(AIcc[4]=-2*model03$logLik+2*teta[4]+(2*teta[4]*(teta[4]+1))/(nrow(data_milkg)-teta[4]-1))

[1] 12224.5

Full WOOD model with class of fixed effect, heterogeneous variance along the time, repeated measures

parf=c()
for (i in 1:length(summary(model03)$ coefficients$fixed))      {
par1=summary(model03)$ coefficients$fixed[[i]]
parf=rbind(parf,par1)
 }

model04=nlme(MY~a*(DIM**b)*exp(-k*DIM),fixed=list(a+b+k~CG-1),random=pdDiag(a+b+k~1),
             control=nlmeControl(maxIter=100),
             data=data_milkg,start=parf,na.action=na.omit,corr=corCAR1(),weights= varPower()) 
parf=c()
for (i in 1:length(summary(model03)$ coefficients$fixed))      {
par1=summary(model03)$ coefficients$fixed[[i]]
parf=rbind(parf,par1)
 }

model04=nlme(MY~a*(DIM**b)*exp(-k*DIM),fixed=list(a+b+k~CG-1),random=pdDiag(a+b+k~1),
             control=nlmeControl(maxIter=100),
             data=data_milkg,start=parf,na.action=na.omit,corr=corCAR1(),weights= varPower()) 

parf=c()
for (i in 1:length(summary(model03)$ coefficients$fixed))      {
par1=summary(model03)$ coefficients$fixed[[i]]
parf=rbind(parf,par1)
 }

model04=nlme(MY~a*(DIM**b)*exp(-k*DIM),fixed=list(a+b+k~CG-1),random=pdDiag(a+b+k~1),
             control=nlmeControl(maxIter=100),
             data=data_milkg,start=parf,na.action=na.omit,corr=corCAR1(),weights= varPower()) 
parf=c()
for (i in 1:length(summary(model03)$ coefficients$fixed))      {
par1=summary(model03)$ coefficients$fixed[[i]]
parf=rbind(parf,par1)
 }

model04=nlme(MY~a*(DIM**b)*exp(-k*DIM),fixed=list(a+b+k~CG-1),random=pdDiag(a+b+k~1),
             control=nlmeControl(maxIter=100),
             data=data_milkg,start=parf,na.action=na.omit,corr=corCAR1(),weights= varPower()) 

parf=c()
for (i in 1:length(summary(model03)$ coefficients$fixed))      {
par1=summary(model03)$ coefficients$fixed[[i]]
parf=rbind(parf,par1)
 }

model04=nlme(MY~a*(DIM**b)*exp(-k*DIM),fixed=list(a+b+k~CG-1),random=pdDiag(a+b+k~1),
             control=nlmeControl(maxIter=100),
             data=data_milkg,start=parf,na.action=na.omit,corr=corCAR1(),weights= varPower()) 
parf=c()
for (i in 1:length(summary(model03)$ coefficients$fixed))      {
par1=summary(model03)$ coefficients$fixed[[i]]
parf=rbind(parf,par1)
 }

model04=nlme(MY~a*(DIM**b)*exp(-k*DIM),fixed=list(a+b+k~CG-1),random=pdDiag(a+b+k~1),
             control=nlmeControl(maxIter=100),
             data=data_milkg,start=parf,na.action=na.omit,corr=corCAR1(),weights= varPower()) 

multi-model selection framework (Burnham and Anderson, 2004)

delta=c()
for(i in 1: length(AIcc)){
delta[i]=AIcc[i]-min(AIcc)
}

wpro=c()
for(i in 1:length(AIcc)){
wpro[i]=exp(-delta[i]/2)
}
sum(wpro[1:length(AIcc)])

[1] 1

wprob=c()
for(i in 1: length(AIcc)){
wprob[i]=exp(-delta[i]/2)/sum(wpro[1:length(AIcc)])
}


ER=c()
for(i in 1: length(AIcc)){
ER[i]=max(wprob)/wprob[i]
}

(quadro.akaike=data.frame(teta,AICc=AIcc,delta,wprob,ER))

Extracting the confidence intervals of the parameters after choosing the best model

intervals(model04)

Prediction MY and TMY

MWoodhat=function(a,b,k,age){
  y=(a*age^b)*exp(-k*age)
  return(y)
}

Total Milk Yield

(TMYCG22=MYWood(a=as.numeric(intervals(model04)$fixed[2,2]),
                 b=as.numeric(intervals(model04)$fixed[8,2]),
                 k=as.numeric(intervals(model04)$fixed[14,2]),
                 MYl=10,MYu=305))

     value   
[1,] 6900.513

(TMYCG32=MYWood(a=as.numeric(intervals(model04)$fixed[5,2]),
               b=as.numeric(intervals(model04)$fixed[11,2]),
               k=as.numeric(intervals(model04)$fixed[17,2]),
               MYl=10,MYu=305))

     value   
[1,] 6739.791

TMYCG22[[1]]-TMYCG32[[1]]

[1] 160.7221

creating interactive graphics with plotly

age=seq(1,305,by=0.25)
#predictions
predc1=MWoodhat(a=as.numeric(intervals(model04)$fixed[2,2]),b=as.numeric(intervals(model04)$fixed[8,2]),k=as.numeric(intervals(model04)$fixed[14,2]),age=age)
predc2=MWoodhat(a=as.numeric(intervals(model04)$fixed[5,2]),b=as.numeric(intervals(model04)$fixed[11,2]),k=as.numeric(intervals(model04)$fixed[17,2]),age=age)

### Multiple Y Axes

#ay <- list(
#  tickfont = list(color = "black"),
#  overlaying = "y",
#  side = "right",
#  title = "Daily Weight Gain(g)"
#)
fig <- plot_ly()
fig <- fig %>% add_lines(x = age, y = predc1, name = "Milk Yield CG22")

fig <- fig %>% add_lines(x = age, y = predc2, name = "Milk Yield CG32")


fig <- fig %>% layout(
  title = "WOOD Milk curve",
  xaxis = list(title="DIM(days)"),
  yaxis = list(title="MY (kg)")
)
fig